function [b,vc,J] = vcqr(bhat,y,x,tau)
% Compute standard errors for IQR when the equation is just identified
% using the method proposed in Powell (1986).  This implementation uses
% a Gaussian kernel and a simple rule of thumb bandwidth.
% [b,vc,J] = vciqr(bhat,y,x,tau)
% Inputs -
%   bhat - estimated coefficients from QR
%   y - dependent variable
%   d - RHS endogenous variable
%   x - covariates (If there are no covariates, pass x = [].)
%   z - instruments
%   tau - quantile index
% Outputs -
%   b - estimated coefficients with standard errors
%   vc - covariance matrix of b
%   J - matrix in asymptotic variance formula (J'SJ) used in process testing


n = size(y,1);	    % Number of observations
x = [x,ones(n,1)];    % Add constant term
k = size(x,2);   % Number of regressors

vc = zeros(k,k);
b = zeros(k,2);
S = (1/n)*x'*x; 

e = y - x*bhat; % Generate residuals

h = 1.364*((2*sqrt(pi))^(-1/5))*sqrt(var(e))*(n^(-1/5)); %Calculate bandwidth using Silverman's rule of thumb
J = (1/(n*h))*((normpdf(e/h)*ones(1,size(x,2))).*x)'*x;

vc = (1/n)*(tau-tau^2)*inv(J')*S*inv(J);
b(:,1) = bhat;
b(:,2) = (sqrt(diag(vc)));

J = inv(J);
